Take_Home Exercise 1: Geospatial Analytics for Social Good: Application of Spatial and Spatio-temporal Point Patterns Analysis to discover the geographical distribution of Armed Conflict in Myanmar

Author

Jiale SO

Published

August 22, 2024

Modified

September 20, 2024

1.0 Introduction

The study of armed conflicts in Myanmar has gained critical importance in understanding the geographical distribution and intensity of violence across different regions. Myanmar’s complex ethnic composition and ongoing civil strife make it a unique case for geospatial analysis. This project aims to apply spatial and spatio-temporal point pattern analysis methods to uncover the patterns of armed conflict between January 2021 and June 2024.

By leveraging conflict data from the Armed Conflict Location & Event Data (ACLED) and geospatial tools, we will focus on visualizing and interpreting conflict density through heat maps, Kernel Density Estimation (KDE), and advanced spatio-temporal analysis.

Our analysis will focus on four types of conflict events:

  1. Battles,
  2. Explosion/Remote violence,
  3. Strategic developments,
  4. Violence against civilians,

with particular attention paid to quarterly patterns in conflict occurrence.

2.0 Setting Up The Environment & Dataset

2.1 Installing the required Packages

Key Packages Used in the Project:

  1. sf: Handles simple features in R, allowing for spatial data manipulation and analysis. It is crucial for reading and managing geospatial data like shapefiles (e.g., Myanmar’s administrative boundaries).

  2. raster: Used for raster-based spatial data manipulation, especially for working with raster maps, such as Kernel Density Estimation (KDE) results.

  3. spatstat: A powerful package for spatial point pattern analysis. It helps to analyze and visualize spatial point data, particularly for identifying clusters or patterns in armed conflict events.

  4. sparr: Builds on spatstat and focuses on performing spatial and spatio-temporal kernel smoothing, which will be crucial for KDE and heatmap creation.

  5. tmap: A thematic mapping package that will allow us to create maps, including KDE visualizations on an OpenStreetMap base.

  6. tidyverse: A collection of data manipulation packages like dplyr, ggplot2, and purrr. It’s essential for data cleaning, manipulation, and visualization tasks.

  7. stpp: Used for spatio-temporal point pattern analysis, crucial for analyzing how conflict events evolve in both space and time.

  8. skimr: A quick and comprehensive tool to provide summaries and descriptive statistics for datasets, helping in the initial exploration.

  9. gganimate: Extends ggplot2 to create animated visualizations. We can use this for animated time-series or evolving conflict maps.

  10. ggplot2: The core plotting package in R, essential for creating visualizations like time series plots and KDE heatmaps.

  11. plotly: Useful for creating interactive visualizations, allowing users to explore spatial data interactively (e.g., hover over points to see conflict details).

  12. pacman: is a package management tool in R designed to streamline the process of loading and installing packages.

pacman::p_load(sf, raster, spatstat, sparr, tmap, tidyverse, stpp, skimr, gganimate, ggplot2, plotly, flexdashboard, shiny, DT,gridExtra, stars)

2.2 Data-set involved in this topic

For this analysis, we use two key datasets:

2.2.1 ACLED Armed Conflict Data

Location & Event Data (ACLED)platform, which maintains an extensive record of conflict events globally. For this specific analysis, we will limit the dataset by filtering based on the following parameters to streamline data preparation and minimize the need for extensive data cleaning:

Data Parameter Filter Category
Date Range From 01/01/2021 to 30/06/2024.
Event Type 1. Battles
2. Violence Against Civilians
3. Explosions/Remote Violence
4. Strategic Developments
Country Myanmar
ACLED Configuration Image

Code to Import ACLED Dataset
ACLEDData <- read_csv("data/raw/aspatial/2021-01-01-2024-06-30-Myanmar.csv")
Rows: 42608 Columns: 35
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (20): event_id_cnty, event_date, disorder_type, event_type, sub_event_ty...
dbl (15): year, time_precision, inter1, inter2, interaction, iso, latitude, ...

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

2.2.1.1 Understanding the data set fields.

Referencing this ACLED Official codebook, this is the dataset that we are working with, not to bore you with the details are mainly interested in the following fields,

  • Event ID: Unique identifier for each conflict event.

  • Event Date: Date of occurrence.

  • Event Type: Type of conflict event (e.g., Battles, Remote Violence).

  • Latitude/Longitude: Coordinates of the event.

  • Fatalities: Number of fatalities resulting from the event.

  • Actors: The groups or individuals involved in the conflict (e.g., state actors, ethnic armed groups).

  • Admin Levels: Administrative region, district, and township where the event took place.

If you’re interested in the data set fields to explore more, here’s the full fields!

ACLED Full Table Fields Summary
Fields name Fields Description Values
event_id_cnty A unique alphanumeric event identifier by number and country acronym. This identifier remains constant even when the event details are updated. E.g. ETH9766
event_date The date on which the event took place. Recorded as Year-Month-Day. E.g. 2023-02-16
year The year in which the event took place. E.g. 2018
time_precision A numeric code between 1 and 3 indicating the level of precision of the date recorded for the event. The higher the number, the lower the precision. 1, 2, or 3; with 1 being the most precise.
disorder_type The disorder category an event belongs to. Political violence, Demonstrations, or Strategic developments.
event_type The type of event; further specifies the nature of the event. E.g. BattlesFor the full list of ACLED event types, see the ACLED Event Types table.
sub_event_type A subcategory of the event type. E.g. Armed clashFor the full list of ACLED sub-event types, see the ACLED Event Types table.
actor1 One of two main actors involved in the event (does not necessarily indicate the aggressor). E.g. Rioters (Papua New Guinea)
assoc_actor_1 Actor(s) involved in the event alongside ‘Actor 1’ or actor designations that further identify ‘Actor 1’. E.g. Labor Group (Spain); Women (Spain)Can have multiple actors separated by a semicolon, or can be blank.
inter1 A numeric code between 0 and 8 indicating the type of ‘Actor 1’ (for more, see the section Actor Names, Types, and ‘Inter’ Codes). 1, 2, 3, 4, 5, 6, 7, or 8.
actor2 One of two main actors involved in the event (does not necessarily indicate the target or victim). E.g. Civilians (Kenya)Can be blank.
assoc_actor_2 Actor(s) involved in the event alongside ‘Actor 2’ or actor designation further identifying ‘Actor 2’. E.g. Labor Group (Spain); Women (Spain)Can have multiple actors separated by a semicolon, or can be blank.
inter2 A numeric code between 0 to 8 indicating the type of ‘Actor 2’ (for more, see the section Actor Names, Types, and ‘Inter’ Codes). 0, 1, 2, 3, 4, 5, 6, 7, or 8.
interaction A two-digit numeric code (combination of ‘Inter 1’ and ‘Inter 2’) indicating the two actor types interacting in the event (for more, see the section Actor Names, Types, and ‘Inter’ Codes). E.g.3, 58
civilian_targeting This column indicates whether the event involved civilian targeting. Either ‘Civilians targeted’ or blank.
iso A unique three-digit numeric code assigned to each country or territory according to ISO 3166. E.g. 231 for Ethiopia
region The region of the world where the event took place. E.g. Eastern Africa
country The country or territory in which the event took place. E.g. Ethiopia
admin1 The largest sub-national administrative region in which the event took place. E.g. Oromia
admin2 The second largest sub-national administrative region in which the event took place. E.g. ArsiCan be blank.
admin3 The third largest sub-national administrative region in which the event took place. E.g. MertiCan be blank.
location The name of the location at which the event took place. E.g. Abomsa
latitude The latitude of the location in four decimal degrees notation (EPSG:32647). E.g. 8.5907
longitude The longitude of the location in four decimal degrees notation (EPSG:32647). E.g. 39.8588
geo_precision A numeric code between 1 and 3 indicating the level of certainty of the location recorded for the event. The higher the number, the lower the precision. 1, 2, or 3; with 1 being the most precise.
source The sources used to record the event. Separated by a semicolon. E.g. Ansar Allah; Yemen Data Project
source_ scale An indication of the geographic closeness of the used sources to the event (for more, see the section Source Scale). E.g. Local partner-National
notes A short description of the event. E.g. On 16 February 2023, OLF-Shane abducted an unidentified number of civilians after stopping a vehicle in an area near Abomsa (Merti, Arsi, Oromia). The abductees were traveling from Adama to Abomsa, Arsi.
fatalities The number of reported fatalities arising from an event. When there are conflicting reports, the most conservative estimate is recorded. E.g. 3No information on fatalities is recorded as 0 reported fatalities.
tags Additional structured information about the event. Separated by a semicolon. E.g. women targeted: politicians; sexual violence
timestamp An automatically generated Unix timestamp that represents the exact date and time an event was uploaded to the ACLED API. E.g. 1676909320

2.2.2 Myanmar Administrative Boundaries (Shapefiles):

Obtained through Geonode Mimu, this shapefile helps us to build the map and set the boundary zone of each district of myanmar. This dataset provides the geographical boundaries of Myanmar’s administrative divisions, from the national level down to the township level. It is essential for mapping conflict events to specific regions.

Note

Myanmar has State, District and Township level, why District Level?

Choosing the district level over the township level for conflict analysis provides a better balance between detail and clarity. The district level allows us to capture broader regional trends without overwhelming the analysis with too many granular data points, as township-level data can be overly detailed. It also improves computational efficiency and makes visualizations clearer, while still offering enough specificity to reveal conflict hotspots. Additionally, population and auxiliary data are more readily available at the district level, making the analysis more consistent and manageable.

Code to Import Shapefile Dataset
M_State_Sf <- st_read(dsn="data/raw/geospatial/stateLevel", layer = "mmr_polbnda_adm1_250k_mimu_1") 
Reading layer `mmr_polbnda_adm1_250k_mimu_1' from data source 
  `C:\Users\jiale\Desktop\IS415\IS415-GAA\Take_Home_Exercises\Take_Home_Exercise_1\data\raw\geospatial\stateLevel' 
  using driver `ESRI Shapefile'
Simple feature collection with 15 features and 6 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 92.1721 ymin: 9.696844 xmax: 101.17 ymax: 28.54554
Geodetic CRS:  WGS 84
M_State_Sf
Simple feature collection with 15 features and 6 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 92.1721 ymin: 9.696844 xmax: 101.17 ymax: 28.54554
Geodetic CRS:  WGS 84
First 10 features:
   OBJECTID          ST ST_PCODE           ST_RG          ST_MMR PCode_V
1         1  Ayeyarwady   MMR017          Region  ဧရာဝတီတိုင်းဒေသကြီး     9.4
2         2        Bago   MMR111          Region    ပဲခူးတိုင်းဒေသကြီး     9.4
3         4        Chin   MMR004           State       ချင်းပြည်နယ်     9.4
4         5      Kachin   MMR001           State       ကချင်ပြည်နယ်     9.4
5         6       Kayah   MMR002           State       ကယားပြည်နယ်     9.4
6         7       Kayin   MMR003           State        ကရင်ပြည်နယ်     9.4
7         8      Magway   MMR009          Region   မကွေးတိုင်းဒေသကြီး     9.4
8         9    Mandalay   MMR010          Region မန္တလေးတိုင်းဒေသကြီး     9.4
9        10         Mon   MMR011           State         မွန်ပြည်နယ်     9.4
10       11 Nay Pyi Taw   MMR018 Union Territory        နေပြည်တော်     9.4
                         geometry
1  MULTIPOLYGON (((95.20798 15...
2  MULTIPOLYGON (((96.17964 19...
3  MULTIPOLYGON (((93.36931 24...
4  MULTIPOLYGON (((97.59674 28...
5  MULTIPOLYGON (((97.1759 19....
6  MULTIPOLYGON (((97.81508 16...
7  MULTIPOLYGON (((94.11699 22...
8  MULTIPOLYGON (((96.14023 23...
9  MULTIPOLYGON (((97.73689 15...
10 MULTIPOLYGON (((96.32013 20...
M_District_Sf <- st_read(dsn="data/raw/geospatial", layer = "mmr_polbnda_adm2_250k_mimu") 
Reading layer `mmr_polbnda_adm2_250k_mimu' from data source 
  `C:\Users\jiale\Desktop\IS415\IS415-GAA\Take_Home_Exercises\Take_Home_Exercise_1\data\raw\geospatial' 
  using driver `ESRI Shapefile'
Simple feature collection with 80 features and 7 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 92.1721 ymin: 9.696844 xmax: 101.17 ymax: 28.54554
Geodetic CRS:  WGS 84
M_District_Sf
Simple feature collection with 80 features and 7 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 92.1721 ymin: 9.696844 xmax: 101.17 ymax: 28.54554
Geodetic CRS:  WGS 84
First 10 features:
   OBJECTID          ST ST_PCODE         DT   DT_PCODE      DT_MMR PCode_V
1         1  Ayeyarwady   MMR017   Hinthada MMR017D002    ဟင်္သာတခရိုင်     9.4
2         2  Ayeyarwady   MMR017    Labutta MMR017D004    လပွတ္တာခရိုင်     9.4
3         3  Ayeyarwady   MMR017     Maubin MMR017D005     မအူပင်ခရိုင်     9.4
4         4  Ayeyarwady   MMR017  Myaungmya MMR017D003 မြောင်းမြခရိုင်     9.4
5         5  Ayeyarwady   MMR017    Pathein MMR017D001      ပုသိမ်ခရိုင်     9.4
6         6  Ayeyarwady   MMR017     Pyapon MMR017D006     ဖျာပုံခရိုင်     9.4
7         7 Bago (East)   MMR007       Bago MMR007D001      ပဲခူးခရိုင်     9.4
8         8 Bago (East)   MMR007    Taungoo MMR007D002    တောင်ငူခရိုင်     9.4
9         9 Bago (West)   MMR008       Pyay MMR008D001      ပြည်ခရိုင်     9.4
10       10 Bago (West)   MMR008 Thayarwady MMR008D002   သာယာဝတီခရိုင်     9.4
                         geometry
1  MULTIPOLYGON (((95.12637 18...
2  MULTIPOLYGON (((95.04462 15...
3  MULTIPOLYGON (((95.38231 17...
4  MULTIPOLYGON (((94.6942 16....
5  MULTIPOLYGON (((94.27572 15...
6  MULTIPOLYGON (((95.20798 15...
7  MULTIPOLYGON (((95.90674 18...
8  MULTIPOLYGON (((96.17964 19...
9  MULTIPOLYGON (((95.70458 19...
10 MULTIPOLYGON (((95.85173 18...

2.2.2.1 Understanding the data set fields.

Field Name Description
OBJECTID This is a unique identifier for each feature in the dataset, typically used to identify individual records or polygons in the shapefile.
ST This represents the State or Region in Myanmar. For example, in your dataset, “Ayeyarwady” refers to a state/region.
ST_PCODE This stands for State Postal Code. It is a standardized code that represents each state or region in Myanmar, such as “MMR017” for Ayeyarwady.
DT This stands for District or Township within the respective state/region. For example, “Hinthada” is a district or township within Ayeyarwady.
DT_PCODE This stands for District/Township Postal Code. It is a standardized postal code for each district or township, such as “MMR017D002” for the Hinthada district/township in Ayeyarwady.
DT_MMR This field could be the District/Township name in Myanmar script, written in the local language. It may be an alternative representation of the “DT” field, showing the name of the district/township in Myanmar’s native language.
PCODE_V This could be a Version of the Postal Code or a verification value used internally in the dataset. In this case, the value is “9.4”, possibly indicating a specific version of postal codes or an accuracy measure.
geometry This column represents the spatial data for each district/township. It contains the geometrical shape (MULTIPOLYGON) defining the boundaries of the state or district/township, with coordinates provided in longitude and latitude.

3.0 Data Pre-processing

To ensure accuracy and usability of the data, several preprocessing steps will be undertaken for the different datasets.

3.1 Myanmar Shapefile

3.1.1 Setting the CRS for the

Since Myanmar uses CRS of 32647 and when we download the map it’s in WGS84, we should change it to 32647 .

st_crs(M_State_Sf)
Coordinate Reference System:
  User input: WGS 84 
  wkt:
GEOGCRS["WGS 84",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["latitude",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["longitude",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    ID["EPSG",4326]]
# Set the CRS for m_sf, assuming the appropriate CRS is WGS 84 (EPSG:32647)
M_State_Sf <- st_transform(M_State_Sf, crs = 32637)
# Verify that the CRS has been correctly set
print(st_crs(M_State_Sf))
Coordinate Reference System:
  User input: EPSG:32637 
  wkt:
PROJCRS["WGS 84 / UTM zone 37N",
    BASEGEOGCRS["WGS 84",
        ENSEMBLE["World Geodetic System 1984 ensemble",
            MEMBER["World Geodetic System 1984 (Transit)"],
            MEMBER["World Geodetic System 1984 (G730)"],
            MEMBER["World Geodetic System 1984 (G873)"],
            MEMBER["World Geodetic System 1984 (G1150)"],
            MEMBER["World Geodetic System 1984 (G1674)"],
            MEMBER["World Geodetic System 1984 (G1762)"],
            MEMBER["World Geodetic System 1984 (G2139)"],
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]],
            ENSEMBLEACCURACY[2.0]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4326]],
    CONVERSION["UTM zone 37N",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",39,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",0.9996,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",500000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Navigation and medium accuracy spatial referencing."],
        AREA["Between 36°E and 42°E, northern hemisphere between equator and 84°N, onshore and offshore. Djibouti. Egypt. Eritrea. Ethiopia. Georgia. Iraq. Jordan. Kenya. Lebanon. Russian Federation. Saudi Arabia. Somalia. Sudan. Syria. Türkiye (Turkey). Ukraine."],
        BBOX[0,36,84,42]],
    ID["EPSG",32637]]
st_crs(M_District_Sf)
Coordinate Reference System:
  User input: WGS 84 
  wkt:
GEOGCRS["WGS 84",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["latitude",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["longitude",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    ID["EPSG",4326]]
# Set the CRS for m_sf, assuming the appropriate CRS is WGS 84 (EPSG:32647)
M_District_Sf <- st_transform(M_District_Sf, crs = 32637)
# Verify that the CRS has been correctly set
print(st_crs(M_District_Sf))
Coordinate Reference System:
  User input: EPSG:32637 
  wkt:
PROJCRS["WGS 84 / UTM zone 37N",
    BASEGEOGCRS["WGS 84",
        ENSEMBLE["World Geodetic System 1984 ensemble",
            MEMBER["World Geodetic System 1984 (Transit)"],
            MEMBER["World Geodetic System 1984 (G730)"],
            MEMBER["World Geodetic System 1984 (G873)"],
            MEMBER["World Geodetic System 1984 (G1150)"],
            MEMBER["World Geodetic System 1984 (G1674)"],
            MEMBER["World Geodetic System 1984 (G1762)"],
            MEMBER["World Geodetic System 1984 (G2139)"],
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]],
            ENSEMBLEACCURACY[2.0]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4326]],
    CONVERSION["UTM zone 37N",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",39,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",0.9996,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",500000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Navigation and medium accuracy spatial referencing."],
        AREA["Between 36°E and 42°E, northern hemisphere between equator and 84°N, onshore and offshore. Djibouti. Egypt. Eritrea. Ethiopia. Georgia. Iraq. Jordan. Kenya. Lebanon. Russian Federation. Saudi Arabia. Somalia. Sudan. Syria. Türkiye (Turkey). Ukraine."],
        BBOX[0,36,84,42]],
    ID["EPSG",32637]]

3.1.2 Renaming and removal of column names

colnames(M_State_Sf) <- c("OBJECTID", "state","state_code","type", "state_myr", "mimi_version", "geometry")

M_State_Sf_Cleansed <- M_State_Sf %>% select(state, type, state_myr ,geometry)
colnames(M_District_Sf) <- c("OBJECTID", "state", "state_code", "district", "district_code", "district_mmr", "mimi_version", "geometry")

M_District_Sf_Cleansed <- M_District_Sf %>% select(state, district, district_mmr, geometry)

3.1.3 Checking for validity of data

When working with spatial data, it’s crucial to ensure that all geometries are valid. Invalid geometries can cause errors in analysis and visualization.

  1. Checking Validity with st_is_valid():

  2. Identifying Invalid Geometries:

  3. Fixing Invalid Geometries with st_make_valid()

#checking if it's valid
M_State_Sf_Validity <- st_is_valid(M_State_Sf_Cleansed)
M_State_Sf_Invalid <- which(!M_State_Sf_Validity)
if (length(M_State_Sf_Invalid) > 0) {
  print("MPZ Invalid!")
  print(M_State_Sf_cleansed[M_State_Sf_Invalid, ])
} else {
  print("it's valid!")
}
[1] "it's valid!"
#checking if it's valid
M_District_Sf_Validity <- st_is_valid(M_District_Sf_Cleansed)
M_District_Sf_Invalid <- which(!M_District_Sf_Validity)
if (length(M_District_Sf_Invalid) > 0) {
  print("MPZ Invalid!")
  print(M_District_Sf_cleansed[M_District_Sf_Invalid, ])
} else {
  print("it's valid!")
}
[1] "it's valid!"

3.1.4 Visualizing the mynamar map

On the top right, you can toggle between the district level and also the state level to understand more about the boundaries of Myanmar.

tmap_mode("plot")
tm_shape(M_State_Sf_Cleansed) +  # Base map (Myanmar boundaries)
  tm_polygons("state",  # Color the base map by state
              palette = "Set3",  # Use Set3 color palette for states
              border.col = "gray",  # Border color for the states
              alpha = 0.5,  # Semi-transparent polygons
              title = "State",  # Legend title for states
              legend.show = TRUE,  # Show legend for state colors
             ) +
  tm_layout(main.title = "States of Myanmar",  # Main map title
            legend.outside = TRUE)  # Position the legend outside the map

There is more than 80 district here, so it only showcases 30 :)

tm_shape(M_District_Sf_Cleansed) +  # Base map (Myanmar boundaries)
  tm_polygons("district",  # Color the base map by district
              palette = "Set3",  # Use Set3 color palette for districts
              border.col = "gray",  # Border color for the districts
              alpha = 0.5,  # Semi-transparent polygons
              title = "District",  # Legend title for districts
              legend.show = TRUE,  # Show legend for district colors
              ) +
  tm_layout(main.title = "Districts of Myanmar",  # Main map title
            legend.outside = TRUE)  # Position the legend outside the map
Warning: Number of levels of the variable "district" is 80, which is larger
than max.categories (which is 30), so levels are combined. Set
tmap_options(max.categories = 80) in the layer function to show all levels.

3.2 ACLED Data

3.2.1 Changing the Column Names

Since Myanmar’s regional hierarchy follows State, District, and Township levels, we will rename the columns accordingly:

  • admin1State

  • admin2District

  • admin3Township

This is important because different countries use different administrative hierarchies. For example, in Singapore, the hierarchy is organized by Region and Subzones. Adjusting these names ensures that our dataset aligns with Myanmar’s specific regional structure for accurate analysis.

ACLEDData_Cleanse <- ACLEDData %>%
  select(event_id_cnty, event_date, year, disorder_type, event_type, actor1, inter1, 
         actor2, inter2, interaction, admin1, admin2, admin3, location, latitude, 
         longitude, fatalities) %>%
  rename(state = admin1, district = admin2, township = admin3) %>%
  mutate(across(where(is.character), ~ replace_na(.x, "NA")),  # Replace NA in character columns with "NA"
         across(where(is.numeric), ~ replace_na(.x, 0)))  # Replace NA in numeric columns with 0

3.2.2 Adding a “Quarter-Year” Column

To facilitate our temporal analysis, we need to add a “quarter-year” column based on the event_date field. This can be done by adjusting the date format to represent the quarter and year, ensuring that each event is categorized by the specific quarter it occurred in (e.g., Q1-2021, Q2-2022). This will allow for easier analysis of conflict trends over time.

# Convert event_date to Date format (if it's not already a date)
ACLEDData_Cleanse$event_date <- as.Date(ACLEDData_Cleanse$event_date, format="%d-%b-%y")  # Adjust the format if needed
# Add a new column that shows the quarter and year
ACLEDData_Cleanse <- ACLEDData_Cleanse %>%
  mutate(quarter_year = paste0("Q", quarter(event_date), "-", year(event_date)))

head(ACLEDData_Cleanse)
# A tibble: 6 × 18
  event_id_cnty event_date  year disorder_type   event_type actor1 inter1 actor2
  <chr>         <date>     <dbl> <chr>           <chr>      <chr>   <dbl> <chr> 
1 MMR64313      2024-06-30  2024 Political viol… Battles    Peopl…      3 Milit…
2 MMR64320      2024-06-30  2024 Political viol… Battles    Peopl…      3 Milit…
3 MMR64321      2024-06-30  2024 Political viol… Battles    Peopl…      3 Milit…
4 MMR64322      2024-06-30  2024 Strategic deve… Strategic… Milit…      1 NA    
5 MMR64323      2024-06-30  2024 Political viol… Battles    PKDF …      3 Milit…
6 MMR64324      2024-06-30  2024 Strategic deve… Strategic… Milit…      1 NA    
# ℹ 10 more variables: inter2 <dbl>, interaction <dbl>, state <chr>,
#   district <chr>, township <chr>, location <chr>, latitude <dbl>,
#   longitude <dbl>, fatalities <dbl>, quarter_year <chr>

3.2.3 Joining ACLED’s Codebook Description

ACLED’s stores their data for the column “interaction” and “inter1” and “inter2” in codes, using their code book, let’s reorganise their data for simplier view, we can reference the code book here to know what each code represent. Map it out as a csv file read it in and change accordingly.

3.2.3.1 Left joining inter1 and inter’s description.

For more details about each inter code read here.

ACLEDActorInterCode <- read_csv("data/raw/aspatial/ActorTypesInterCode.csv")
Rows: 8 Columns: 2
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (1): Description
dbl (1): code

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
ACLEDData_Cleanse <- ACLEDData_Cleanse %>%
  left_join(ACLEDActorInterCode, by = c("inter1" = "code")) %>%
  rename(inter1_description = Description)
# Left join again for inter2
ACLEDData_Cleanse <- ACLEDData_Cleanse %>%
  left_join(ACLEDActorInterCode, by = c("inter2" = "code")) %>%
  rename(inter2_description = Description)

datatable(head(ACLEDData_Cleanse, 5), options = list(pageLength = 5, autoWidth = TRUE))

3.2.3.2 Left joining interaction description.

For more details about each interaction code read here.

ACLEDInteractionCode <- read_csv("data/raw/aspatial/AcledInteractionCodes.csv")
Rows: 44 Columns: 2
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (1): Description
dbl (1): code

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
ACLEDData_Cleanse <- ACLEDData_Cleanse %>%
  left_join(ACLEDInteractionCode, by = c("interaction" = "code")) %>%
  rename(interaction_description = Description)

head(ACLEDData_Cleanse)
# A tibble: 6 × 21
  event_id_cnty event_date  year disorder_type   event_type actor1 inter1 actor2
  <chr>         <date>     <dbl> <chr>           <chr>      <chr>   <dbl> <chr> 
1 MMR64313      2024-06-30  2024 Political viol… Battles    Peopl…      3 Milit…
2 MMR64320      2024-06-30  2024 Political viol… Battles    Peopl…      3 Milit…
3 MMR64321      2024-06-30  2024 Political viol… Battles    Peopl…      3 Milit…
4 MMR64322      2024-06-30  2024 Strategic deve… Strategic… Milit…      1 NA    
5 MMR64323      2024-06-30  2024 Political viol… Battles    PKDF …      3 Milit…
6 MMR64324      2024-06-30  2024 Strategic deve… Strategic… Milit…      1 NA    
# ℹ 13 more variables: inter2 <dbl>, interaction <dbl>, state <chr>,
#   district <chr>, township <chr>, location <chr>, latitude <dbl>,
#   longitude <dbl>, fatalities <dbl>, quarter_year <chr>,
#   inter1_description <chr>, inter2_description <chr>,
#   interaction_description <chr>

3.2.3 Making it a SF Object and reverse geolocate state and district

Since ACLED provides longitude and latitude data, I prefer to reverse geolocate the points to match Myanmar’s official state and district boundaries. We are uncertain how ACLED assigns these regions, so to ensure consistency, we remove the original state and district columns from the ACLED data and replace them with the geolocated values.

Steps:

  1. Convert ACLED Data to an SF Object: Using longitude and latitude coordinates, transform the ACLED dataset into a spatial object. Remember taht we have to set CRS 32647 here as well.

  2. Perform a Spatial Join: Match the points from ACLED with the corresponding state and district boundaries from the m_sf shapefile, selecting only those columns.

  3. Remove Original Columns: After the spatial join, drop the original state, district, and township columns from the ACLED dataset.

  4. Rename the Joined Columns: Rename the newly joined state.y and district.y to state and district, effectively replacing the original columns with the reverse-geolocated values.

# Step 1: Convert ACLEDDataCleanse to an sf object using longitude and latitude
ACLEDData_Cleanse_Sf <- st_as_sf(ACLEDData_Cleanse, coords = c("longitude", "latitude"), crs = 4326, remove = FALSE)

# Step 2: Transform ACLEDDataCleanse_sf to the same CRS as m_sf (EPSG: 32637)
ACLEDData_Cleanse_Sf <- st_transform(ACLEDData_Cleanse_Sf, crs = 32637)

# Step 3: Perform a spatial join, selecting only the state and district from m_sf
reverse_geolocated_state <- st_join(ACLEDData_Cleanse_Sf, M_State_Sf_Cleansed[, c("state")], join = st_intersects)

# Step 2: Spatial join to add 'district' from M_District_Sf_Cleansed
reverse_geolocated_district <- st_join(reverse_geolocated_state, M_District_Sf_Cleansed[, c("district")], join = st_intersects)

# Step 3: Remove original 'state', 'district', and 'township' columns from ACLEDData_Cleanse (if they exist)
# This step removes the original columns, and then renames the newly joined columns
ACLEDData_Cleanse_Sf <- reverse_geolocated_district %>%
  select(-contains("state.x"), -contains("district.x"), -contains("township")) %>%  # Remove old state, district, township columns
  rename(state = state.y, district = district.y)  # Rename newly joined columns

ACLEDData_Cleanse <- st_drop_geometry(ACLEDData_Cleanse_Sf)

# View the updated data
print(ACLEDData_Cleanse_Sf)
Simple feature collection with 42608 features and 20 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 6629596 ymin: 2135072 xmax: 8594728 ymax: 5015434
Projected CRS: WGS 84 / UTM zone 37N
# A tibble: 42,608 × 21
   event_id_cnty event_date  year disorder_type  event_type actor1 inter1 actor2
   <chr>         <date>     <dbl> <chr>          <chr>      <chr>   <dbl> <chr> 
 1 MMR64313      2024-06-30  2024 Political vio… Battles    Peopl…      3 Milit…
 2 MMR64320      2024-06-30  2024 Political vio… Battles    Peopl…      3 Milit…
 3 MMR64321      2024-06-30  2024 Political vio… Battles    Peopl…      3 Milit…
 4 MMR64322      2024-06-30  2024 Strategic dev… Strategic… Milit…      1 NA    
 5 MMR64323      2024-06-30  2024 Political vio… Battles    PKDF …      3 Milit…
 6 MMR64324      2024-06-30  2024 Strategic dev… Strategic… Milit…      1 NA    
 7 MMR64325      2024-06-30  2024 Political vio… Battles    Milit…      1 PSLF/…
 8 MMR64326      2024-06-30  2024 Political vio… Battles    PSLF/…      2 Milit…
 9 MMR64328      2024-06-30  2024 Political vio… Battles    Milit…      1 PSLF/…
10 MMR64330      2024-06-30  2024 Political vio… Battles    Milit…      1 PSLF/…
# ℹ 42,598 more rows
# ℹ 13 more variables: inter2 <dbl>, interaction <dbl>, location <chr>,
#   latitude <dbl>, longitude <dbl>, fatalities <dbl>, quarter_year <chr>,
#   inter1_description <chr>, inter2_description <chr>,
#   interaction_description <chr>, geometry <POINT [m]>, state <chr>,
#   district <chr>

3.2.5 Visualizing it by Event Type

# Set tmap mode to plot for static maps
tmap_mode("plot")
tmap mode set to plotting
bbox <- st_bbox(M_State_Sf_Cleansed)  # Assuming you are using sf object

# Create the tmap object with the base map and event markers
tm_shape(M_State_Sf_Cleansed) + 
    tm_polygons("state", alpha = 0.3, border.col = "gray") +  # Base map with higher transparency
  
  # Add event markers (bubbles) with size based on fatalities
  tm_shape(ACLEDData_Cleanse_Sf) + 
  tm_bubbles(size = "fatalities",  # Marker size based on fatalities
             col = "event_type",  # Color markers by event type
             palette = "Set1",  # Use Set1 color palette for event types
             border.col = "black",  # Border color for bubbles
             border.alpha = 0.5,  # Semi-transparent border
             title.size = "Number of Fatalities",  # Title for bubble size legend
             title.col = "Event Types") +  # Title for event type legend
  
  # Layout settings for the map, including title and legend positioning
  tm_layout(main.title = "Myanmar's State Conflicts by Fatalities",  # Main map title
            main.title.size = 1.5,  # Title font size
            legend.outside = TRUE,  # Position legend outside the map
            legend.outside.size = 0.5,  # Adjust size of outside legend
            legend.position = c("left", "top"),  # Position for the event type legend
            legend.title.size = 1.2,  # Size of the legend title
            legend.text.size = 1,  # Size of the legend text
            legend.bg.color = "white",  # Background color for the legend
            legend.bg.alpha = 0.5,  # Transparency for the legend background
            inner.margins = c(0.05, 0.05, 0.05, 0.05)) 

4.0 Exploratory Data analysis

With the data cleansed, let’s conduct a high-level analysis of our dataset to determine some key statistics.

The central question we’re exploring is:

“Which district in Myanmar has the highest number of conflicts that may be concerning to civilians?”

In section 3.2.5, we recognized that Myanmar has many states and numerous conflicts during this period. Therefore, we will aggregate the data at the state level to prioritize states based on the following key figures:

  • Civilian Involvement: We focus on events where either Actor 1 or Actor 2 is a civilian, as these are likely to raise humanitarian concerns.

  • Total Conflicts: This measures the overall number of conflict events in each state, helping us identify which states experienced the most conflict.

  • Conflict Density: This calculates the number of conflicts per square kilometer in each state, allowing us to understand the intensity of conflicts relative to the size of the state.

  • Total Fatalities: This sums up the total number of fatalities within each state, highlighting areas with the highest death toll, which are likely experiencing the most severe impacts.

  • Fatalities Density: This measures the fatalities per square kilometer, giving a sense of how deadly the conflicts are in each state in relation to its area.

This analysis will allow us to focus on the states most affected by conflicts that pose the greatest risk to civilians.

4.1 Aggregation of Data for exploratory purposes

Step 1: Filter Data for Civilians

In this first step, the data (ACLEDData_Cleanse) is filtered so that only rows where either “Civilians” are involved in the conflict (inter1_description or inter2_description) are kept. Additionally, rows where the state information is missing (NA) are removed because these conflicts can’t be assigned to a specific geographical area.

filtered_data <- ACLEDData_Cleanse %>%
  filter((inter1_description == "Civilians" | inter2_description == "Civilians") & !is.na(state))

Step 2: Calculate State Area in Square Kilometers

Next, the code calculates the area for each state from the M_State_Sf_Cleansed spatial dataset. The st_area function retrieves the area in square meters, and dividing by 1 million (1e6) converts it into square kilometers. This information will be used later to calculate conflict and fatality density.

state_area_km2 <- st_area(M_State_Sf_Cleansed) / 1e6
state_area_df <- data.frame(state = M_State_Sf_Cleansed$state, area_km2 = as.numeric(state_area_km2))

Step 3: Aggregate Total Conflicts and Fatalities by State

The next step is to group the filtered data by state and calculate two key metrics: the total number of conflicts and the total number of fatalities for each state. This is done using group_by and summarise.

agg_data_by_state <- filtered_data %>%
  group_by(state) %>%
  summarise(
    total_conflicts = n(), 
    total_fatalities = sum(fatalities, na.rm = TRUE), 
    .groups = "drop"
  )

Step 4: Convert to Non-Spatial Data Frame

Here, the spatial data is converted into a regular data frame by dropping the geometry information. This allows easier manipulation of the non-spatial data for further processing.

filtered_data_df <- st_drop_geometry(filtered_data)

Step 5: Prepare Data Frame for Yearly Data

To store yearly conflict and fatality data for each state, an empty data frame is initialized by selecting distinct states from the filtered data

final_result <- filtered_data_df %>% select(state) %>% distinct()

Step 6: Loop Through Each Year and Calculate Yearly Conflicts/Fatalities

In this step, the code iterates over each unique year in the dataset. For each year, it calculates the total conflicts and fatalities for each state. The results are merged back into the final_result data frame, with columns named according to the year (e.g., conflicts_2021, fatalities_2021).

unique_years <- unique(filtered_data_df$year)

for (yr in unique_years) {
  yearly_data <- filtered_data_df %>%
    filter(year == yr) %>%
    group_by(state) %>%
    summarise(
      conflicts = n(), 
      fatalities = sum(fatalities, na.rm = TRUE)
    ) %>%
    rename_at(vars(conflicts, fatalities), ~ paste0(., "_", yr))  # Rename with year
  
  final_result <- left_join(final_result, yearly_data, by = "state")
}

Step 7: Replace NAs with 0

Any missing data for a state-year combination (e.g., a state with no conflicts in a particular year) is replaced with 0 to avoid leaving gaps in the analysis.

final_result <- final_result %>%
  mutate(across(everything(), ~ ifelse(is.na(.), 0, .)))

Step 8: Calculate Conflict and Fatality Density

Now, the code calculates the density of conflicts and fatalities per square kilometer for each state. The results from Step 3 (agg_data_by_state) are joined with the state area data from Step 2. The density is simply the number of conflicts or fatalities divided by the area in square kilometers.

agg_data_with_density <- agg_data_by_state %>%
  left_join(state_area_df, by = "state") %>%
  mutate(
    conflict_density = total_conflicts / area_km2,
    fatality_density = total_fatalities / area_km2
  )

Step 9: Merge Yearly Data and Density with Spatial Data

Finally, the spatial data (M_State_Sf_Cleansed) is merged with the yearly conflict/fatality data (final_result) and the density data (agg_data_with_density). Additional calculations are made to find the fatality-per-conflict ratio which are then added to the final spatial dataset.

finalized_map <- M_State_Sf_Cleansed %>%
  left_join(final_result, by = "state") %>%
  left_join(agg_data_with_density, by = "state")

finalized_map <- finalized_map %>%
  mutate(conflict_fatality_ratio = ifelse(total_fatalities == 0, NA, total_fatalities / total_conflicts))

finalized_map <- finalized_map %>%
  mutate(
    fatality_per_conflict = ifelse(total_conflicts == 0, NA, total_fatalities / total_conflicts),
  )

finalized_map <- finalized_map %>%
  mutate(
    rank_total_conflicts = rank(-total_conflicts, ties.method = "min"),  # Rank total_conflicts (largest to smallest)
    rank_total_fatalities = rank(-total_fatalities, ties.method = "min"),  # Rank total_fatalities (largest to smallest)
    rank_conflict_density = rank(-conflict_density, ties.method = "min"),  # Rank conflict_density (largest to smallest)
    rank_fatality_density = rank(-fatality_density, ties.method = "min"),  # Rank fatality_density (largest to smallest)
    rank_fatality_per_conflict = rank(-fatality_per_conflict, ties.method = "min")  # Rank fatality_per_conflict (largest to smallest)
  )

Final Output: Display the Data

The finalized_map contains all the processed information, including conflict/fatality counts, densities, and state geometry, which can now be visualized or further analyzed.

head(finalized_map)
Simple feature collection with 6 features and 23 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 6589238 ymin: 2946007 xmax: 8116377 ymax: 5139186
Projected CRS: WGS 84 / UTM zone 37N
       state   type      state_myr conflicts_2024 fatalities_2024
1 Ayeyarwady Region ဧရာဝတီတိုင်းဒေသကြီး             64               9
2       Bago Region   ပဲခူးတိုင်းဒေသကြီး            194             128
3       Chin  State      ချင်းပြည်နယ်             72              52
4     Kachin  State      ကချင်ပြည်နယ်            160              95
5      Kayah  State      ကယားပြည်နယ်             21              38
6      Kayin  State       ကရင်ပြည်နယ်             52              59
  conflicts_2023 fatalities_2023 conflicts_2022 fatalities_2022 conflicts_2021
1             83               7            177              29            233
2            339             220            172              61            191
3            101              76            204              80            160
4            316             133            226             140            221
5            113              74            149              87            114
6            165              83            165              86            116
  fatalities_2021 total_conflicts total_fatalities  area_km2 conflict_density
1              34             557               79  91740.37      0.006071482
2              87             896              496 106296.56      0.008429247
3              32             537              240  83773.78      0.006410120
4              71             923              439 214569.60      0.004301635
5             108             397              307  33443.39      0.011870806
6              57             498              285  92162.25      0.005403514
  fatality_density                       geometry conflict_fatality_ratio
1     0.0008611258 MULTIPOLYGON (((7515288 297...               0.1418312
2     0.0046661904 MULTIPOLYGON (((7376294 368...               0.5535714
3     0.0028648583 MULTIPOLYGON (((6594137 416...               0.4469274
4     0.0020459562 MULTIPOLYGON (((6707242 513...               0.4756230
5     0.0091796915 MULTIPOLYGON (((7484860 384...               0.7732997
6     0.0030923725 MULTIPOLYGON (((7888506 331...               0.5722892
  fatality_per_conflict rank_total_conflicts rank_total_fatalities
1             0.1418312                   11                    14
2             0.5535714                    9                     6
3             0.4469274                   12                    13
4             0.4756230                    8                     7
5             0.7732997                   14                    10
6             0.5722892                   13                    12
  rank_conflict_density rank_fatality_density rank_fatality_per_conflict
1                    11                    15                         15
2                     8                     8                          4
3                    10                    10                          9
4                    14                    12                          6
5                     7                     4                          1
6                    12                     9                          3

4.2 Visualizing the aggregated data.

For each column that we created earlier, we can now showcase it in a choropleth map to highlight the states with the highest values in each category. These maps will help visually identify which state has the highest value for each column. You can navigate through the tabset to explore the different metrics and easily view the data on a map!

# Create the tmap object for total civilian conflicts
tm_total_conflicts <- tm_basemap(server = c("Esri.WorldGrayCanvas", "OpenStreetMap", "Esri.WorldTopoMap")) +
  # Layer for total civilian conflicts
  tm_shape(finalized_map) +
  tm_polygons("total_conflicts", 
              title = "Total Civilian Conflicts",
              palette = "Reds", 
              border.col = "black", 
              alpha = 0.8,  # Transparency for polygons
              border.alpha = 0.5,  # Semi-transparent border
              id = "state") +  # Set state as identifier for the polygons

  # Layout settings matching your style
  tm_layout(
    main.title = "Myanmar's State by Total Civilian Related Conflicts",  # Main map title
    main.title.size = 1.5,  # Title font size
    legend.outside = TRUE,  # Position legend outside the map
    legend.outside.size = 0.5,  # Adjust size of outside legend
    legend.position = c("left", "top"),  # Position for the event type legend
    legend.title.size = 1.2,  # Size of the legend title
    legend.text.size = 1,  # Size of the legend text
    legend.bg.color = "white",  # Background color for the legend
    legend.bg.alpha = 0.5,  # Transparency for the legend background
    inner.margins = c(0.05, 0.05, 0.05, 0.05)  # Set inner margins for better spacing
  )
tm_total_conflicts

# Create the tmap object for total civilian conflicts
tm_total_fatalities <- tm_basemap(server = c("Esri.WorldGrayCanvas", "OpenStreetMap", "Esri.WorldTopoMap")) +
  # Layer for total civilian conflicts
  tm_shape(finalized_map) +
  tm_polygons("total_fatalities", 
              title = "Total Fatalities by Civilian Related Conflicts",
              palette = "Blues", 
              border.col = "black", 
              alpha = 0.8,  # Transparency for polygons
              border.alpha = 0.5,  # Semi-transparent border
              id = "state") +  # Set state as identifier for the polygons

  # Layout settings matching your style
  tm_layout(
    main.title = "Myanmar's State by Total Fatatlies For Civilian Related Conflicts",  # Main map title
    main.title.size = 1.0,  # Title font size
    legend.outside = TRUE,  # Position legend outside the map
    legend.outside.size = 0.5,  # Adjust size of outside legend
    legend.position = c("left", "top"),  # Position for the event type legend
    legend.title.size = 1.2,  # Size of the legend title
    legend.text.size = 1,  # Size of the legend text
    legend.bg.color = "white",  # Background color for the legend
    legend.bg.alpha = 0.5,  # Transparency for the legend background
    inner.margins = c(0.05, 0.05, 0.05, 0.05)  # Set inner margins for better spacing
  )

tm_total_fatalities

tm_conflict_density <- tm_basemap(server = c("Esri.WorldGrayCanvas", "OpenStreetMap", "Esri.WorldTopoMap")) +
  # Layer for total civilian conflicts
  tm_shape(finalized_map) +
  tm_polygons("conflict_density", 
              title = "Myanmar State by Civilian Related Conflicts Density",
              palette = "Blues", 
              border.col = "black", 
              alpha = 0.8,  # Transparency for polygons
              border.alpha = 0.5,  # Semi-transparent border
              id = "state") +  # Set state as identifier for the polygons

  # Layout settings matching your style
  tm_layout(
    main.title = "Myanmar State by Civilian Related Conflicts Density",  # Main map title
    main.title.size = 1.0,  # Title font size
    legend.outside = TRUE,  # Position legend outside the map
    legend.outside.size = 0.5,  # Adjust size of outside legend
    legend.position = c("left", "top"),  # Position for the event type legend
    legend.title.size = 1.2,  # Size of the legend title
    legend.text.size = 1,  # Size of the legend text
    legend.bg.color = "white",  # Background color for the legend
    legend.bg.alpha = 0.5,  # Transparency for the legend background
    inner.margins = c(0.05, 0.05, 0.05, 0.05)  # Set inner margins for better spacing
  )

tm_conflict_density

tm_fatalities_conflict <- tm_basemap(server = c("Esri.WorldGrayCanvas", "OpenStreetMap", "Esri.WorldTopoMap")) +
  # Layer for total civilian conflicts
  tm_shape(finalized_map) +
  tm_polygons("fatality_density", 
              title = "Myanmar State by Fatalies For Civilians Related Conflicts Density Per Km^2",
              palette = "Blues", 
              border.col = "black", 
              alpha = 0.8,  # Transparency for polygons
              border.alpha = 0.5,  # Semi-transparent border
              id = "state") +  # Set state as identifier for the polygons

  # Layout settings matching your style
  tm_layout(
    main.title = "Myanmar State by Fatalies For Civilians Related Conflicts Density Per Km^2",  # Main map title
    main.title.size = 1.0,  # Title font size
    legend.outside = TRUE,  # Position legend outside the map
    legend.outside.size = 0.5,  # Adjust size of outside legend
    legend.position = c("left", "top"),  # Position for the event type legend
    legend.title.size = 1.0,  # Size of the legend title
    legend.text.size = 0.8,  # Size of the legend text
    legend.bg.color = "white",  # Background color for the legend
    legend.bg.alpha = 0.5,  # Transparency for the legend background
    inner.margins = c(0.05, 0.05, 0.05, 0.05)  # Set inner margins for better spacing
  )

tm_fatalities_conflict

tm_fatalities_density <- tm_basemap(server = c("Esri.WorldGrayCanvas", "OpenStreetMap", "Esri.WorldTopoMap")) +
  # Layer for total civilian conflicts
  tm_shape(finalized_map) +
  tm_polygons("fatality_per_conflict", 
              title = "Myanmar State by Fatalies for Civilian Related Conflicts",
              palette = "Blues", 
              border.col = "black", 
              alpha = 0.8,  # Transparency for polygons
              border.alpha = 0.5,  # Semi-transparent border
              id = "state") +  # Set state as identifier for the polygons

  # Layout settings matching your style
  tm_layout(
    main.title = "Myanmar State by Fatalies for Civilian Related Conflicts",  # Main map title
    main.title.size = 1.0,  # Title font size
    legend.outside = TRUE,  # Position legend outside the map
    legend.outside.size = 0.5,  # Adjust size of outside legend
    legend.position = c("left", "top"),  # Position for the event type legend
    legend.title.size = 1.0,  # Size of the legend title
    legend.text.size = 0.8,  # Size of the legend text
    legend.bg.color = "white",  # Background color for the legend
    legend.bg.alpha = 0.5,  # Transparency for the legend background
    inner.margins = c(0.05, 0.05, 0.05, 0.05)  # Set inner margins for better spacing
  )

tm_fatalities_density

4.3 Analysis on the aggregated data to find top 3 state to analyse

4.3.1 How to choose which States?

With the maps above, everyone can choose what they would like to analyze, whether it’s total civilian conflicts or fatalities by state. However, I would like to focus on the central question:

Which state has the highest amount of conflict and fatality density?

By using density, I’ve aimed to assess fatalities in relation to conflicts and ensure a fair comparison by adjusting for density.

Note

Keep in mind, this is not the most precise analysis, as other factors—such as ethnicity, population demographics, and gender—could also be considered. However, for the purpose of this study, I have chosen to focus on conflict and fatality density.

Conflict Density: This is the number of conflicts per square kilometer in a state.

  1. Conflict Density=Total ConflictsArea of State (km²) = Conflict Density=Area of State (km²)Total Conflicts​

    This tells you how many conflicts occur per unit of area (1 km²) in each state.

  2. Fatality Density: This is the number of fatalities per square kilometer.

    Fatality Density=Total FatalitiesArea of State (km²) = Fatality Density=Area of State (km²)Total Fatalities​

4.3.2 Viewing the top 3 states across

Using the map below, you can interactively explore the Conflict and Fatality Density. We observe that the top three states are as follows, and here’s how they compare across the data set.

You can use the layer toggle to switch between fatalities and conflicts, or click on the map layers to view detailed, aggregated data for each state.

# Step 1: Set tmap mode to view for interactive maps
tmap_mode("view")
tmap mode set to interactive viewing
# Step 2: Define popup variables to show all relevant columns when clicking on a state
popup_variables <- c( 
                     "Total Civilian Conflicts" = "total_conflicts",
                     "Fatalities" = "total_fatalities",
                     "Civilian Conflict Density" = "conflict_density",
                     "Fatality Density" = "fatality_density",
                     "Fatality Per Civilian Conflict" = "fatality_per_conflict",
                     "Ranked For Total Conflict By State (Out of 15)" = "rank_total_conflicts",
                     "Ranked For Total Fatalities By State (Out of 15)" = "rank_total_fatalities",
                     "Ranked For Conflict Density By State (Out of 15)" = "rank_conflict_density",
                     "Ranked For Fatalities Density By State (Out of 15)" = "rank_fatality_density",
                     "Ranked For Fatalies/Conflict By State (Out of 15)" = "rank_fatality_per_conflict"
                     )

# Step 3: Create the map with 6 layers for toggling
tm <- tm_basemap(server = c("Esri.WorldGrayCanvas", "OpenStreetMap", "Esri.WorldTopoMap")) +
  tm_shape(finalized_map) +
  tm_polygons("conflict_density", 
              title = "Civilian Conflict Density",
              palette = "Oranges", 
              border.col = "black",
              popup.vars = popup_variables,
              id = "state", 
              group = "Civilian Conflict Density") +
  tm_shape(finalized_map) +
  tm_polygons("fatality_density", 
              title = "Fatality Density",
              palette = "Purples", 
              border.col = "black",
              popup.vars = popup_variables,
              id = "state", 
              group = "Fatality Density") +
  tm_layout(
    legend.outside = TRUE,
    legend.outside.size = 0.5,  # Adjust the size of the legend
    legend.position = c("left", "top")  # Position the legend
  )
# Step 5: Display the interactive map
tm
legend.postion is used for plot mode. Use view.legend.position in tm_view to set the legend position in view mode.
tmap_mode("plot")
tmap mode set to plotting
finalized_data_df <- as.data.frame(st_drop_geometry(finalized_map))
# Sort the data by conflict_density and fatality_density in descending order
finalized_data_df <- finalized_data_df %>%
  arrange(desc(conflict_density), desc(fatality_density))

# Format the densities to show "per km²"
finalized_data_df$conflict_density_label <- paste0(round(finalized_data_df$conflict_density, 2), " per km²")
finalized_data_df$fatality_density_label <- paste0(round(finalized_data_df$fatality_density, 2), " per km²")

# Create bar chart for Conflict Density (sorted by conflict_density)
conflict_density_plot <- ggplot(finalized_data_df, aes(x = reorder(state, -conflict_density), y = conflict_density)) +
  geom_bar(stat = "identity", fill = "orange") +
  geom_text(aes(label = conflict_density_label), vjust = 1.8, size = 1.5, color = "black", position = position_stack(vjust = 0.5)) +  # Place text at bottom of bar
  ggtitle("Top States by Conflict Density (per km²)") +
  xlab("State") +
  ylab("Conflict Density (per km²)") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

# Create bar chart for Fatality Density (sorted by fatality_density)
fatality_density_plot <- ggplot(finalized_data_df, aes(x = reorder(state, -fatality_density), y = fatality_density)) +
  geom_bar(stat = "identity", fill = "purple") +
  geom_text(aes(label = fatality_density_label), vjust = 1.8, size = 1.5, color = "black", position = position_stack(vjust = 0.5)) +  # Place text at bottom of bar
  ggtitle("Top States by Fatality Density (per km²)") +
  xlab("State") +
  ylab("Fatality Density (per km²)") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

# Display both charts side by side
grid.arrange(conflict_density_plot, fatality_density_plot, ncol = 2)

4.3.3 Top 3 Chosen state for further analysis

4.3.3.1 Yangon

4.3.3.2 Mandalay

4.3.3.3 Sagaling.

4.4 Final Round of Cleansing ACLED Data

Now that we have finalized the data and narrowed down the states that we are going to analyse, we will filter out for the states data.

Yangon_ACLED_Data <- filtered_data %>% 
  filter(state %in% c("Yangon"))

class(Yangon_ACLED_Data)
[1] "tbl_df"     "tbl"        "data.frame"
Yangon_ACLED_Data_Sf <- st_as_sf(Yangon_ACLED_Data, coords = c("longitude", "latitude"), crs = 4326, remove = FALSE)

Yangon_ACLED_Data_Sf <- st_transform(Yangon_ACLED_Data_Sf, crs = 32637)


class(Yangon_ACLED_Data_Sf)
[1] "sf"         "tbl_df"     "tbl"        "data.frame"
Yangon_ACLED_Validity <- st_is_valid(Yangon_ACLED_Data_Sf)
Yangon_invalid <- which(!Yangon_ACLED_Validity)
if (length(Yangon_invalid) > 0) {
  print("Yangon is Invalid!")
  print(Yangon_ACLED_Data_Sf[Yangon_invalid, ])
} else {
  print("Yangon_ACLED_Data_Sf is valid!")
}
[1] "Yangon_ACLED_Data_Sf is valid!"
Yangon_ACLED_SFDF <- as_Spatial(Yangon_ACLED_Data_Sf)
Yangon_ACLED_SFDF
class       : SpatialPointsDataFrame 
features    : 1489 
extent      : 7470244, 7680669, 3162567, 3345642  (xmin, xmax, ymin, ymax)
crs         : +proj=utm +zone=37 +datum=WGS84 +units=m +no_defs 
variables   : 20
names       : event_id_cnty, event_date, year,          disorder_type,                 event_type,                       actor1, inter1,                actor2, inter2, interaction,      location, latitude, longitude, fatalities, quarter_year, ... 
min values  :      MMR10949,      18633, 2021,     Political violence, Explosions/Remote violence,           5 Brothers Younger,      1, Civilians (Australia),      0,          17,     Ah Hpyauk,  16.4389,   95.6934,          0,      Q1-2021, ... 
max values  :      MMR64089,      19895, 2024, Strategic developments, Violence against civilians, YUG: Yangon Urban Guerrillas,      7,                    NA,      7,          70, Zee Hpyu Kone,  17.5612,   96.7448,         13,      Q4-2023, ... 
Mandalay_ACLED_Data <- filtered_data %>% 
  filter(state %in% c("Mandalay"))

class(Mandalay_ACLED_Data)
[1] "tbl_df"     "tbl"        "data.frame"
Mandalay_ACLED_Data_Sf <- st_as_sf(Mandalay_ACLED_Data, coords = c("longitude", "latitude"), crs = 4326, remove = FALSE)

class(Mandalay_ACLED_Data_Sf)
[1] "sf"         "tbl_df"     "tbl"        "data.frame"
Mandalay_ACLED_Validity <- st_is_valid(Mandalay_ACLED_Data_Sf)
Mandalay_invalid <- which(!Mandalay_ACLED_Validity)
if (length(Mandalay_invalid) > 0) {
  print("Mandalay is Invalid!")
  print(Mandalay_ACLED_Data_Sf[Mandalay_invalid, ])
} else {
  print("Mandalay_ACLED_Data_Sf is valid!")
}
[1] "Mandalay_ACLED_Data_Sf is valid!"
Mandalay_ACLED_SFDF <- as_Spatial(Mandalay_ACLED_Data_Sf)
Mandalay_ACLED_SFDF
class       : SpatialPointsDataFrame 
features    : 2021 
extent      : 94.849, 96.6148, 20.4319, 23.667  (xmin, xmax, ymin, ymax)
crs         : +proj=longlat +datum=WGS84 +no_defs 
variables   : 20
names       : event_id_cnty, event_date, year,          disorder_type,                 event_type,                                   actor1, inter1,              actor2, inter2, interaction,         location, latitude, longitude, fatalities, quarter_year, ... 
min values  :      MMR10960,      18640, 2021,     Political violence, Explosions/Remote violence, 21KPG: 21 Guerrilla Force - Kyaukpadaung,      1, Civilians (Myanmar),      0,          17, 4 Maing Kan Thar,  20.4319,    94.849,          0,      Q1-2021, ... 
max values  :      MMR64318,      19904, 2024, Strategic developments, Violence against civilians,          Zero Guerrilla Force - Myingyan,      7,                  NA,      7,          70,       Zin Chaung,   23.667,   96.6148,         10,      Q4-2023, ... 
Sagaing_ACLED_Data <- filtered_data %>% 
  filter(state %in% c("Sagaing"))

class(Sagaing_ACLED_Data)
[1] "tbl_df"     "tbl"        "data.frame"
Sagaing_ACLED_Data_Sf <- st_as_sf(Sagaing_ACLED_Data, coords = c("longitude", "latitude"), crs = 4326, remove = FALSE)

class(Sagaing_ACLED_Data_Sf)
[1] "sf"         "tbl_df"     "tbl"        "data.frame"
Sagaing_ACLED_Validity <- st_is_valid(Sagaing_ACLED_Data_Sf)
Sagaing_invalid <- which(!Sagaing_ACLED_Validity)
if (length(Sagaing_invalid) > 0) {
  print("Sagaing is Invalid!")
  print(Sagaing_ACLED_Data_Sf[Sagaing_invalid, ])
} else {
  print("Sagaing_ACLED_Data_Sf is valid!")
}
[1] "Sagaing_ACLED_Data_Sf is valid!"
Sagaing_ACLED_SFDF <- as_Spatial(Sagaing_ACLED_Data_Sf)
Sagaing_ACLED_SFDF
class       : SpatialPointsDataFrame 
features    : 5346 
extent      : 93.9575, 96.6034, 21.605, 26.3006  (xmin, xmax, ymin, ymax)
crs         : +proj=longlat +datum=WGS84 +no_defs 
variables   : 20
names       : event_id_cnty, event_date, year,          disorder_type,                 event_type,                                      actor1, inter1,            actor2, inter2, interaction,               location, latitude, longitude, fatalities, quarter_year, ... 
min values  :      MMR10952,      18640, 2021,     Political violence, Explosions/Remote violence, ABSDF: All Burma Students' Democratic Front,      1, Civilians (China),      0,          17,               55 Maing,   21.605,   93.9575,          0,      Q1-2021, ... 
max values  :      MMR65891,      19898, 2024, Strategic developments, Violence against civilians,             Zero Guerrilla Force - Myingyan,      7,                NA,      7,          70, Zin Ka Le (Zin Ka Lin),  26.3006,   96.6034,        175,      Q4-2023, ... 

Viola now we have the final state that that en

5.0 Spatial Data Wrangling

Now that we have chosen our 3 states to analyse we have to prepare the data to analyse them within

5.1 Setting up the Owin Window for States.

In this part of the code, we are setting up a spatial window for Myanmar, which defines the boundaries within which the spatial point pattern analysis will take place. This is done using the owin function from the spatstat package, which converts the polygonal boundary of Myanmar into a format that can be used for further spatial analysis.

Yangon_Sf <- M_State_Sf %>% 
  filter(state %in% c("Yangon"))
# Step 1: Extract the individual polygons from the multipolygon
yangon_polygons <- st_cast(Yangon_Sf$geometry, "POLYGON")
yangon_polygons_filtered <- yangon_polygons[-c(1, 2)]
yangon_multipolygon_filtered <- st_combine(yangon_polygons_filtered)
# Add the filtered multipolygon back to the Yangon_Sf object
Yangon_Sf$geometry <- yangon_multipolygon_filtered
Yangon <- as_Spatial(Yangon_Sf)
Yangon_Owin <- as.owin(Yangon_Sf)
plot(Yangon_Owin)

summary(Yangon_Owin)
Window: polygonal boundary
6 separate polygons (1 hole)
                 vertices         area relative.area
polygon 1            3832  2.79639e+10      1.00e+00
polygon 2 (hole)        3 -5.16800e+02     -1.85e-08
polygon 3              14  6.05403e+05      2.16e-05
polygon 4              11  4.56071e+05      1.63e-05
polygon 5              20  5.15327e+06      1.84e-04
polygon 6              37  7.35365e+06      2.63e-04
enclosing rectangle: [7458583, 7698916] x [3150830, 3397886] units
                     (240300 x 247100 units)
Window area = 27977500000 square units
Fraction of frame area: 0.471
Mandalay_Sf <- M_State_Sf %>% 
  filter(state %in% c("Mandalay"))
Mandalay <- as_Spatial(Mandalay_Sf)
Mandalay_Owin <- as.owin(Mandalay_Sf)
plot(Mandalay_Owin)

summary(Mandalay_Owin)
Window: polygonal boundary
single connected closed polygon with 5914 vertices
enclosing rectangle: [6990555, 7369406] x [3760343, 4336524] units
                     (378900 x 576200 units)
Window area = 79199300000 square units
Fraction of frame area: 0.363
Sagaing_Sf <- M_State_Sf %>% 
  filter(state %in% c("Sagaing"))
Sagaing <- as_Spatial(Sagaing_Sf)
Sagaing_Owin <- as.owin(Sagaing_Sf)
plot(Sagaing_Owin)

summary(Sagaing_Owin)
Window: polygonal boundary
single connected closed polygon with 5882 vertices
enclosing rectangle: [6600929, 7143229] x [3910303, 4901317] units
                     (542300 x 991000 units)
Window area = 2.21637e+11 square units
Fraction of frame area: 0.412

5.2 Setting up the ACLED Spatial Class

5.2.1 Converting Spatial Class Into Generic PPP Format

# Yangon_ACLED_Data
# Yangon_ACLED_Data_Sf
# Yangon_ACLED_SFDF 

Yangon_ACLED_coords <- st_coordinates(Yangon_ACLED_Data_Sf)

# Define the window using the bounding box
Yangon_ACLED_bbox <- st_bbox(Yangon_ACLED_Data_Sf)
Yangon_ACLED_window <- owin(xrange = Yangon_ACLED_bbox[c("xmin", "xmax")], yrange = Yangon_ACLED_bbox[c("ymin", "ymax")])

# Create the ppp object
Yangon_ACLED_ppp <- ppp(x = Yangon_ACLED_coords[, 1], y = Yangon_ACLED_coords[, 2], window = Yangon_ACLED_window)
Warning: data contain duplicated points
# Check the ppp object
summary(Yangon_ACLED_ppp)
Planar point pattern:  1489 points
Average intensity 3.865152e-08 points per square unit

*Pattern contains duplicated points*

Coordinates are given to 9 decimal places

Window: rectangle = [7470244, 7680669] x [3162567, 3345642] units
                    (210400 x 183100 units)
Window area = 38523700000 square units
plot(Yangon_ACLED_ppp)

#Mandalay_ACLED_Data
# Mandalay_ACLED_Data_Sf
# Mandalay_ACLED_SFDF 
Mandalay_ACLED_coords <- st_coordinates(Mandalay_ACLED_Data_Sf)

# Define the window using the bounding box
Mandalay_ACLED_bbox <- st_bbox(Mandalay_ACLED_Data_Sf)
Mandalay_ACLED_window <- owin(xrange = Mandalay_ACLED_bbox[c("xmin", "xmax")], yrange = Mandalay_ACLED_bbox[c("ymin", "ymax")])

# Create the ppp object
Mandalay_ACLED_ppp <- ppp(x = Mandalay_ACLED_coords[, 1], y = Mandalay_ACLED_coords[, 2], window = Mandalay_ACLED_window)
Warning: data contain duplicated points
# Check the ppp object
summary(Mandalay_ACLED_ppp)
Planar point pattern:  2021 points
Average intensity 353.7831 points per square unit

*Pattern contains duplicated points*

Coordinates are given to 4 decimal places

Window: rectangle = [94.849, 96.6148] x [20.4319, 23.667] units
                    (1.766 x 3.235 units)
Window area = 5.71254 square units
plot(Mandalay_ACLED_ppp)

Sagaing_ACLED_coords <- st_coordinates(Sagaing_ACLED_Data_Sf)

# Define the window using the bounding box
Sagaing_ACLED_bbox <- st_bbox(Sagaing_ACLED_Data_Sf)
Sagaing_ACLED_window <- owin(xrange = Sagaing_ACLED_bbox[c("xmin", "xmax")], yrange = Sagaing_ACLED_bbox[c("ymin", "ymax")])

# Create the ppp object
Sagaing_ACLED_ppp <- ppp(x = Sagaing_ACLED_coords[, 1], y = Sagaing_ACLED_coords[, 2], window = Sagaing_ACLED_window)
Warning: data contain duplicated points
# Check the ppp object
summary(Sagaing_ACLED_ppp)
Planar point pattern:  5346 points
Average intensity 430.2932 points per square unit

*Pattern contains duplicated points*

Coordinates are given to 4 decimal places

Window: rectangle = [93.9575, 96.6034] x [21.605, 26.3006] units
                    (2.646 x 4.696 units)
Window area = 12.4241 square units
plot(Sagaing_ACLED_ppp)

5.2.2 Dealing with duplicate conflicts

5.2.2.1 Check For Duplicate

any(duplicated(Yangon_ACLED_ppp))
[1] TRUE
sum(multiplicity(Yangon_ACLED_ppp) > 1)
[1] 1432
any(duplicated(Yangon_ACLED_ppp))
[1] TRUE
sum(multiplicity(Yangon_ACLED_ppp) > 1)
[1] 1432

5.2.2.2 Unique Marking them

Yangon_marks <- rep(1, npoints(Yangon_ACLED_ppp))
Yangon_marks[duplicated(Yangon_ACLED_ppp)] <- 2

# Create a new ppp object with marks attached to each point
marked_Yangon_ACLED_ppp <- ppp(x = Yangon_ACLED_coords[, 1], y = Yangon_ACLED_coords[, 2], window = Yangon_ACLED_window, marks = Yangon_marks)
Warning: data contain duplicated points
# Check the ppp object with unique marks
summary(marked_Yangon_ACLED_ppp)
Marked planar point pattern:  1489 points
Average intensity 3.865152e-08 points per square unit

*Pattern contains duplicated points*

Coordinates are given to 9 decimal places

marks are numeric, of type 'double'
Summary:
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  1.000   2.000   2.000   1.916   2.000   2.000 

Window: rectangle = [7470244, 7680669] x [3162567, 3345642] units
                    (210400 x 183100 units)
Window area = 38523700000 square units
plot(marked_Yangon_ACLED_ppp, main = "Yangon Conflicts with Unique Marks", cols = c("black", "red"))

Note

Why use mark and not jitter or delete for these duplicates?

When dealing with conflict data, such as ACLED, uniquely marking events is essential to preserve both the spatial and temporal accuracy of key occurrences. Events often take place at the same location but during different time periods, and jittering the data would distort their true geographic positions, while deleting them would erase important historical information. By marking events as unique, we ensure that the dataset retains its integrity, allowing for the analysis of recurring conflicts without compromising precision. This approach ensures that both spatial and temporal patterns of conflict hotspots are maintained, which is crucial for understanding the dynamics of the events over time.

5.3.1 Combining Point Events and State’s Owin Object

marked_Yangon_ACLED_ppp[Yangon_Owin]
Marked planar point pattern: 1489 points
marks are numeric, of storage type  'double'
window: polygonal boundary
enclosing rectangle: [7458583, 7698916] x [3150830, 3397886] units
plot(marked_Yangon_ACLED_ppp[Yangon_Owin])

5.4 Combining them both

# Step 1: Define the bounding box and window for Yangon
Yangon_ACLED_bbox <- st_bbox(Yangon_Sf)
Yangon_ACLED_window <- owin(
  xrange = c(Yangon_ACLED_bbox["xmin"], Yangon_ACLED_bbox["xmax"]),
  yrange = c(Yangon_ACLED_bbox["ymin"], Yangon_ACLED_bbox["ymax"])
)

# Extract unique quarters
quarters <- unique(Yangon_ACLED_Data_Sf$quarter_year)

# Initialize an empty list to store ppp objects
plot_list <- list()

# Step 2: Loop over each quarter, jitter duplicates, and store ppp objects
for (quarter in quarters) {
  
  # Filter the data for the current quarter
  quarter_data <- Yangon_ACLED_Data_Sf %>% 
    filter(quarter_year == quarter)
  
  # Extract coordinates for this quarter
  quarter_coords <- st_coordinates(quarter_data)
  
  # Create a ppp object for this quarter
  quarter_ppp <- ppp(
    x = quarter_coords[, 1],
    y = quarter_coords[, 2],
    window = Yangon_ACLED_window
  )
  
  # Step 3: Jitter duplicates
  is_duplicate <- duplicated(quarter_ppp)
  jittered_x <- quarter_coords[, 1] + ifelse(is_duplicate, runif(npoints(quarter_ppp), -0.1, 0.1), 0)
  jittered_y <- quarter_coords[, 2] + ifelse(is_duplicate, runif(npoints(quarter_ppp), -0.1, 0.1), 0)
  
  # Create a jittered ppp object with the new coordinates
  jittered_conflict_ppp <- ppp(
    x = jittered_x,
    y = jittered_y,
    window = Yangon_ACLED_window
  )
  
  # Step 4: Subset the jittered points within Yangon window (Yangon owins)
  quarter_conflict_ppp <- jittered_conflict_ppp[Yangon_Owin]
  
  # Step 5: Store the ppp object for later plotting
  plot_list[[quarter]] <- quarter_conflict_ppp
}
Warning: data contain duplicated points
Warning: data contain duplicated points
Warning: data contain duplicated points
Warning: data contain duplicated points
Warning: data contain duplicated points
Warning: data contain duplicated points
Warning: data contain duplicated points
Warning: data contain duplicated points
Warning: data contain duplicated points
Warning: data contain duplicated points
Warning: data contain duplicated points
Warning: data contain duplicated points
Warning: data contain duplicated points
Warning: data contain duplicated points
# Step 6: Plot all quarter data using spatstat's plot function
par(mfrow = c(4, 4))  # Adjust rows and columns for multiple plots

for (quarter in names(plot_list)) {
  plot(plot_list[[quarter]], main = paste("Yangon Conflict Events -", quarter))
}

Yangon_ACLED_bbox <- st_bbox(Yangon_Sf)
Yangon_ACLED_window <- owin(
  xrange = c(Yangon_ACLED_bbox["xmin"], Yangon_ACLED_bbox["xmax"]),
  yrange = c(Yangon_ACLED_bbox["ymin"], Yangon_ACLED_bbox["ymax"])
)

# Extract unique quarters
quarters <- unique(Yangon_ACLED_Data_Sf$quarter_year)

# Initialize an empty dictionary-like list to store ppp objects by quarter
yangon_quarters_conflict <- list()
# Step 2: Loop over each quarter, jitter duplicates, and store the ppp object in the list
for (quarter in quarters) {
  # Filter the data for the current quarter
  quarter_data <- Yangon_ACLED_Data_Sf %>%
    filter(quarter_year == quarter)
  # Extract coordinates for this quarter
  quarter_coords <- st_coordinates(quarter_data)
  # Create a ppp object for this quarter
  quarter_ppp <- ppp(
    x = quarter_coords[, 1],
    y = quarter_coords[, 2],
    window = Yangon_ACLED_window
  )
  # Step 3: Jitter duplicates
  is_duplicate <- duplicated(quarter_ppp)
  jittered_x <- quarter_coords[, 1] + ifelse(is_duplicate, runif(npoints(quarter_ppp), -0.1, 0.1), 0)
  jittered_y <- quarter_coords[, 2] + ifelse(is_duplicate, runif(npoints(quarter_ppp), -0.1, 0.1), 0)
  
  # Create a new ppp object with jittered points
  jittered_ppp <- ppp(
    x = jittered_x,
    y = jittered_y,
    window = Yangon_ACLED_window
  )
  
  # Step 4: Store the jittered ppp object in the dictionary (list) with the quarter as the key
  yangon_quarters_conflict[[quarter]] <- jittered_ppp
}
Warning: data contain duplicated points
Warning: data contain duplicated points
Warning: data contain duplicated points
Warning: data contain duplicated points
Warning: data contain duplicated points
Warning: data contain duplicated points
Warning: data contain duplicated points
Warning: data contain duplicated points
Warning: data contain duplicated points
Warning: data contain duplicated points
Warning: data contain duplicated points
Warning: data contain duplicated points
Warning: data contain duplicated points
Warning: data contain duplicated points
kde_plot_list <- list()

# Step 2: Loop over each quarter, jitter duplicates, and store KDE objects
for (quarter in quarters) {
  
  # Filter the data for the current quarter
  quarter_data <- Yangon_ACLED_Data_Sf %>% 
    filter(quarter_year == quarter)
  
  # Extract coordinates for this quarter
  quarter_coords <- st_coordinates(quarter_data)
  
  # Create a ppp object for this quarter
  quarter_ppp <- ppp(
    x = quarter_coords[, 1],
    y = quarter_coords[, 2],
    window = Yangon_ACLED_window
  )
  
  # Step 3: Jitter duplicates
  is_duplicate <- duplicated(quarter_ppp)
  jittered_x <- quarter_coords[, 1] + ifelse(is_duplicate, runif(npoints(quarter_ppp), -0.1, 0.1), 0)
  jittered_y <- quarter_coords[, 2] + ifelse(is_duplicate, runif(npoints(quarter_ppp), -0.1, 0.1), 0)
  
  # Create a jittered ppp object with the new coordinates
  jittered_conflict_ppp <- ppp(
    x = jittered_x,
    y = jittered_y,
    window = Yangon_ACLED_window
  )
  
  # Step 4: Subset the jittered points within Yangon window (Yangon owins)
  quarter_conflict_ppp <- jittered_conflict_ppp[Yangon_Owin]
  
  # Step 5: Rescale the ppp object to kilometers for KDE
  jittered_ppp_km <- rescale(quarter_conflict_ppp, 1000, "km")
  
  # Step 6: Calculate KDE using Gaussian kernel and bw.scott for bandwidth
  kde_quarter <- density(jittered_ppp_km, sigma = bw.scott, edge = TRUE, kernel = "gaussian")
  
  # Step 7: Store the KDE object for later plotting
  kde_plot_list[[quarter]] <- kde_quarter
}
Warning: data contain duplicated points
Warning: data contain duplicated points
Warning: data contain duplicated points
Warning: data contain duplicated points
Warning: data contain duplicated points
Warning: data contain duplicated points
Warning: data contain duplicated points
Warning: data contain duplicated points
Warning: data contain duplicated points
Warning: data contain duplicated points
Warning: data contain duplicated points
Warning: data contain duplicated points
Warning: data contain duplicated points
Warning: data contain duplicated points
# Step 8: Plot all quarter KDE data using spatstat's plot function
par(mfrow = c(4, 4))  # Adjust rows and columns for multiple plots

for (quarter in names(kde_plot_list)) {
  plot(kde_plot_list[[quarter]], main = paste("KDE for Yangon -", quarter))
}

5.5 Summary of Data sets

We have created a total of five datasets, each serving a distinct purpose in our spatial analysis:

  1. ACLEDData_SF: This is the cleaned and formatted ACLED dataset, stored as an sf (simple feature) object. It contains the key event data, including attributes like location, event type, and time period, and will be used for further spatial analysis.

  2. M_SF: This dataset represents the Myanmar boundary as an sf object. It will be used as the base layer in the analysis, defining the spatial extent and providing a reference map for the region.

  3. Aggregated Level: This dataset is used for exploratory analysis and contains aggregated data for different event types, quarters, or other relevant categories. It allows for a broader overview before diving into more detailed spatial point pattern analyses.

  4. ACLED_ppp_corrected: This is the point pattern dataset created from the ACLED data. It includes the corrected coordinates (with jittering applied) and is stored as a ppp object (planar point pattern) for use in spatstat. This object will be used for more in-depth spatial analyses, such as density estimation or clustering.

  5. myanmar_owin: This is the spatial window object representing the boundary of Myanmar. It defines the geographic limits for our spatial point pattern analysis and ensures that all event data is analyzed within this boundary.

5.6 getting the bandwidth

6.0 Deriving the Quarterly Kernel Density Estimation Layers

6.1 Setting the band widths

# Load necessary libraries
library(spatstat)
library(dplyr)
library(sf)

# Step 1: Define the function to create jittered ppp objects for conflict data
process_quarter_conflicts <- function(region_sf, region_window, region_name, data_sf) {
  
  # Extract unique quarters
  quarters <- unique(data_sf$quarter_year)
  
  # Initialize empty lists to store ppp and KDE objects
  region_quarters_conflict <- list()
  kde_plot_list <- list()
  
  # Step 2: Loop over each quarter and process conflict data
  for (quarter in quarters) {
    
    # Filter the data for the current quarter
    quarter_data <- data_sf %>%
      filter(quarter_year == quarter)
    
    # Extract coordinates for this quarter
    quarter_coords <- st_coordinates(quarter_data)
    
    # Create a ppp object for this quarter
    quarter_ppp <- ppp(
      x = quarter_coords[, 1],
      y = quarter_coords[, 2],
      window = region_window
    )
    
    # Step 3: Remove rejected points that fall outside the window
    valid_points <- inside.owin(quarter_ppp$x, quarter_ppp$y, region_window)
    quarter_ppp <- quarter_ppp[valid_points]
    
    # Step 4: Jitter duplicates
    is_duplicate <- duplicated(quarter_ppp)
    jittered_x <- quarter_ppp$x + ifelse(is_duplicate, runif(npoints(quarter_ppp), -0.1, 0.1), 0)
    jittered_y <- quarter_ppp$y + ifelse(is_duplicate, runif(npoints(quarter_ppp), -0.1, 0.1), 0)
    
    # Create a jittered ppp object with the new coordinates
    jittered_ppp <- ppp(
      x = jittered_x,
      y = jittered_y,
      window = region_window
    )
    
    # Step 5: Store the jittered ppp object for later use
    region_quarters_conflict[[quarter]] <- jittered_ppp
    
    # Step 6: Rescale to kilometers and calculate KDE
    jittered_ppp_km <- rescale(jittered_ppp, 1000, "km")
    kde_quarter <- density(jittered_ppp_km, sigma = 100, edge = TRUE, kernel = "gaussian")
    
    # Step 7: Store the KDE object
    kde_plot_list[[quarter]] <- kde_quarter
  }
  
  # Return the processed ppp and KDE objects
  return(list(
    "ppp_list" = region_quarters_conflict,
    "kde_list" = kde_plot_list
  ))
}

# Step 2: Define Yangon-specific variables
Yangon_ACLED_bbox <- st_bbox(Yangon_ACLED_Data_Sf)
Yangon_ACLED_window <- owin(
  xrange = c(Yangon_ACLED_bbox["xmin"], Yangon_ACLED_bbox["xmax"]),
  yrange = c(Yangon_ACLED_bbox["ymin"], Yangon_ACLED_bbox["ymax"])
)

# Step 3: Process the conflict data for Yangon
results <- process_quarter_conflicts(
  region_sf = Yangon_Sf,
  region_window = Yangon_Owin,
  region_name = "Yangon",
  data_sf = Yangon_ACLED_Data_Sf
)
Warning: data contain duplicated points
Warning: data contain duplicated points
Warning: data contain duplicated points
Warning: data contain duplicated points
Warning: data contain duplicated points
Warning: data contain duplicated points
Warning: data contain duplicated points
Warning: data contain duplicated points
Warning: data contain duplicated points
Warning: data contain duplicated points
Warning: data contain duplicated points
Warning: data contain duplicated points
Warning: data contain duplicated points
Warning: data contain duplicated points
# Step 4: Plot the results
# Plotting the jittered conflict points for each quarter
par(mfrow = c(4, 4))  # Adjust rows and columns for multiple plots

for (quarter in names(results$ppp_list)) {
  plot(results$ppp_list[[quarter]], main = paste("Yangon Conflict Events -", quarter))
}

# Plotting the Kernel Density Estimate (KDE) for each quarter
par(mfrow = c(4, 4))  # Adjust rows and columns for multiple plots

for (quarter in names(results$kde_list)) {
  plot(results$kde_list[[quarter]], main = paste("KDE for Yangon -", quarter))
}

Mandalay_ACLED_bbox <- st_bbox(Mandalay_ACLED_Data_Sf)
Mandalay_ACLED_window <- owin(
  xrange = c(Mandalay_ACLED_bbox["xmin"], Mandalay_ACLED_bbox["xmax"]),
  yrange = c(Mandalay_ACLED_bbox["ymin"], Mandalay_ACLED_bbox["ymax"])
)

# Step 3: Process the conflict data for Mandalay
results <- process_quarter_conflicts(
  region_sf = Mandalay_Sf,
  region_window = Mandalay_Owin,
  region_name = "Mandalay",
  data_sf = Mandalay_ACLED_Data_Sf
)
Warning: 139 points were rejected as lying outside the specified window
Warning: 123 points were rejected as lying outside the specified window
Warning: 152 points were rejected as lying outside the specified window
Warning: 157 points were rejected as lying outside the specified window
Warning: 120 points were rejected as lying outside the specified window
Warning: 139 points were rejected as lying outside the specified window
Warning: 124 points were rejected as lying outside the specified window
Warning: 113 points were rejected as lying outside the specified window
Warning: 149 points were rejected as lying outside the specified window
Warning: 149 points were rejected as lying outside the specified window
Warning: 234 points were rejected as lying outside the specified window
Warning: 161 points were rejected as lying outside the specified window
Warning: 210 points were rejected as lying outside the specified window
Warning: 51 points were rejected as lying outside the specified window
# Step 4: Plot the results
# Plotting the jittered conflict points for each quarter
par(mfrow = c(4, 4))  # Adjust rows and columns for multiple plots

for (quarter in names(results$ppp_list)) {
  plot(results$ppp_list[[quarter]], main = paste("Mandalay Conflict Events -", quarter))
}

# Plotting the Kernel Density Estimate (KDE) for each quarter
par(mfrow = c(4, 4))  # Adjust rows and columns for multiple plots

for (quarter in names(results$kde_list)) {
  plot(results$kde_list[[quarter]], main = paste("KDE for Mandalay -", quarter))
}

6.0 Perform 2nd-Order Spatial Point Pattern Analysis:

7.0 Derive Quarterly Spatio-Temporal KDE Layers:

8.0 Perform 2nd-Order Spatio-Temporal Point Patterns Analysis:

7.0 Results and Interpretation

8.0 Conclusion

9.0 References